% revision from f_run20110927_mixed_sub3items_1figs.m 
% Objective: revised codes around ln 312, so the spatial position is
% modulated by LIP_sal


% function [n t]=f_run0715_mixed_sub3items(Img, fig_opt)
function [n t]=f_run20110403_mixed_sub3items_1fig(Img, fig_opt, strategy, confidence_thresh, noise_sigma, resp_sigma)
flag_subitizing=0;
strategy='mixed';
load('sal_output.mat');
fig_opt='';
%Img=EyeField(:,:);
confidence_thresh= .9;
 noise_sigma=.05;
 resp_sigma=.05;
 title_text='';
 sal_hit_thresh=.6;
 sal_terminate_thresh=.4;
%fig_opt=[];
close all;
init;
field_size=max(size(Img));
NLIP_RANGE=30;
EyeField=Img;
center=[1+floor(field_size/2), 1+floor(field_size/2)];
center=int16(center);
LIP_Eyepos=zeros(field_size);
LIP_Eyepos(1+floor(field_size/2), 1+floor(field_size/2))=0.6;
total_move=100;
% Layers DLPFC_L2 and DLPFC_L1 compose a recurrent network
% which memorize the location of the items.
DLPFC_L2=zeros(field_size);
DLPFC_L1=zeros(field_size);
DLPFC_L2_fr=zeros(field_size);
DLPFC_L1_fr=zeros(field_size);
DLPFC_inh_mask=zeros(field_size);
DLPFC_inh=zeros(field_size);
DLPFC_Low_THRESHOLD=1.4;
DLPFC_L2_thresh = .1;
LIP_sal=zeros(field_size);
LIP_sal_fr=zeros(field_size);
NLIP=zeros(field_size);
XNLIP=zeros(field_size);
NLIP_inh = zeros(field_size);

% DLPFC Global excitation signal send 
DLPFC_HitFromIT = zeros(field_size);
eta_DLPFC_L2=5;
V2=zeros(field_size);
% confidence_thresh=.5;
% noise_sigma=.01;
% DLPFC_NHi_confidence_sigma;

% removied to fit scene integration
% % NLIP confidence mask
% NLIP_resp_mask=zeros(NLIP_RANGE);
% NLIP_c_idx = int16((NLIP_RANGE+1)/2);
% for i=1:NLIP_RANGE
%     for j=1:NLIP_RANGE
%         if ((i-NLIP_c_idx)^2+(j-NLIP_c_idx)^2)<=((NLIP_RANGE-1)/2)^2
%             NLIP_resp_mask(i,j)=1;
%         end
%     end
% end

% removied to fit scene integration
% % PFC_N(prefrontal) rep of numbers 0-6
size_PFC_N=7;
num_count=0:size_PFC_N-1;
% PFC_NR=zeros(1,size_PFC_N);
% PFC_NS=zeros(1,size_PFC_N);
% PFC_NC=zeros(1,size_PFC_N);
% PFC_NR_fr=zeros(1,size_PFC_N);
% PFC_NS_fr=zeros(1,size_PFC_N);
% PFC_NC_fr=zeros(1,size_PFC_N);
% PFC_NShift=zeros(1,size_PFC_N); % neuron with threshold
% shift_thresh=1;
% shift_amp = 50;
% PFC_NC(1,1)=0.6;
% PFC_NC(1,1)=1;

% for diffenetial equation
dt=.08;
duration_per_move=1.2;
pause_time=.1;

    clim_dlpfc=[0 2];
    clim_vin=[.8 1];
    clim_lip_ctr=[0 2];
    pos=center;mv=[0 0];
    stop_flag=0;  % stop the simulation when stop_flag == 1
    act_record=[];

    %% Initially assume V2<-EyeField
    sim_clk=1;
    if ~isequal(fig_opt, 'fig_off')
        figure; title(title_text); pause(3);
        close all;
        h_fig=figure('position',[0 0 600 760], 'color', 'white'); title(title_text);
%                         subplot(5,2,2);bar(num_count,PFC_NS);ylabel('DLPFC-N2');
%                         subplot(5,2,4);bar(num_count,PFC_NR);ylabel('DLPFC-N1');
                        subplot(5,2,6);mesh(NLIP);zlabel('LIP-N');
                        subplot(5,2,9);imagesc(1-EyeField, clim_vin);ylabel('Eye Field');colormap('hot');pause(pause_time);
                        subplot(5,2,5); mesh(LIP_sal_fr);ylabel('LIP/FEF Saliency');pause(pause_time);
                        subplot(5,2,7);mesh(LIP_Eyepos);ylabel('LIP Eye Position');  pause(pause_time);  %imagesc(LIP_Eyepos, clim_lip_ctr);
                        subplot(5,2,3); mesh(DLPFC_L1);ylabel('DLPFC-L1');
                        subplot(5,2,1); (DLPFC_L2);ylabel('DLPFC-L2');
%                  h_fig=figure('position',[-1300 -200 600 1200]); %for portait display
    end
    

    %%
    for i=1:total_move
%   dedicate in tracking one object4, avoid distracted by others
%   return inhibition only when target is recognized

% removied to fit scene integration
%          if i<=3 || strcmp(strategy, 'sub_ohnly')
%              flag_subitizing=1; 
%          end

%% Setting subitizing terminating condition         
%        % Revised sub termination conds; 
%        % subitizing in two cond: 
%           1. i==4 && strategy=='mixed'; 1a. condidence>=confidence, end prog;
%           1b. start counging
%           2. strategy=='sub_only' 2a. condidence>=confidence; 2b. i>8

% removied to fit scene integration
%         [v PFC_N_idx]=max(PFC_NR);
%         confidence = v/sum(PFC_NR);
%         if (strcmp(strategy, 'mixed')&&i==5)
%             flag_subitizing=0;
%             if PFC_N_idx~=1 && confidence>=confidence_thresh
%                  n=min(6,PFC_N_idx-1);
%                  break;
%             else             
%                  PFC_NS=zeros(1,size_PFC_N);
%                  PFC_NR=zeros(1,size_PFC_N);
%                  PFC_NR(1,1)=.7;             
%             end
%         elseif (strcmp(strategy, 'sub_only'))
%             
%             if i<=8
%                 if PFC_N_idx~=1 && confidence>=confidence_thresh
%                     n=min(6,PFC_N_idx-1);
%                     flag_subitizing=0;
%                     break;                    
%                 end
%             else
%                  n=min(6,PFC_N_idx-1);
%                  flag_subitizing=0;
%                  break;                
%             end
%         end
        SumLatInh=0;   %DLPFC_N layer
         

         %^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
         % decide saccade direction vector 'mv'
         % 1. saccade to max of NLIP in subitizing stage
         % 2. saccade to max of LIP_sal in counting stage
         pos=Extract2d_max(LIP_Eyepos);
         
% removied to fit scene integration         
%          % Subitizing
%          if flag_subitizing==1
%             if i==1
%                 tmp=-(LIP_sal_fr(center(1)-1:center(1)+1,center(2)-1:center(2)+1)); % extract min neighboring position 
%                 tmp(2,2)=-1;
%                 target=Extract2d_max(tmp);
%                 try
%                     mv=target-center;
%                 catch
%                     mv=[1 1];
%                 end
%             else
%                 mv=-mv;
%             end
%             mv=mv./abs(mv);
%          % counting
%          else 


            % Find the current max activation in the eye field
            % * has to be implemented as winner take all later.
%             pos=Extract2d_max(LIP_Eyepos);
            mv=Extract_NN(LIP_sal_fr, center, sal_terminate_thresh);
            % catch the no target exception
            try
                mv=mv-center;
            catch
                mv=[1 1];
            end
            mv=mv./abs(mv);
            mv=int16(mv);
%         end
        %^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
        
        %+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
        % at the end of saccade, move eye center away from target to
        % avoid
        % keeping activate LIP neurons
        if stop_flag==1
            %mv=[1 1]; % eye move out the target area
            mv=int16(mv);
        end
        % Update eyefield (changed to update by image and current position)
        % to get current position, pos must plus mv
        EyeField=update_eyefield(Img,pos+mv);
        if ~isequal(fig_opt, 'fig_off')
            subplot(5,2,9);imagesc(1-EyeField, clim_vin);ylabel('Eye Field');colormap('hot');pause(pause_time);
        end
        %+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
        
        % Update V2
        V2=EyeField;
        
        % % compute V2 stimulus to LIP_sal
        % % % this is temporary, since it should either directly saccade to
        % % % destination or get LIP_center info.
        V2_LIP_offset=zeros(field_size);
        V2_LIP_offset(pos(1)+mv(1),pos(2)+mv(2))=1;
        V2_out = conv2(V2, V2_LIP_offset, 'same');
        
% removied to fit scene integration
% %       once hit object(object not in subitizing area), send hit signal to LIP and DLPFC_N
%         if (LIP_sal_fr(center(1), center(2))>=sal_hit_thresh&&NLIP_inh(pos(1),pos(2))==0&&flag_subitizing==0)
%             DLPFC_HitFromIT = ones(field_size);
%             PFC_N_shift_flag=1; % ready to incriment number
% %           EyeField(center(1), center(2))=0;
%         else
%             DLPFC_HitFromIT = zeros(field_size);
%             % update DLPFC_N filter, no input filter
%             PFC_NMask=[0 1 0];
%             PFC_N_shift_flag=0;
%         end
        
        
        % Makes mask that mimics the impulse response of LIP_Eyepos
        w_mask=zeros(3); w_mask(mv(1)+2, mv(2)+2)= 1;
        try
        SubMapLIP = LIP_Eyepos(-1+pos(1):pos(1)+1, -1+pos(2):pos(2)+1);  %the impulse response is very local
        SubMapLIP_fr=mp2fr(SubMapLIP);
        catch
            msg='illegal crop of LIP_Eyepos!'
        end
        
        % Find current peak in DLPFC_N
        % [v, PFC_N_idx] = max(PFC_NHi);

        
        % Dynamic Update
        for j=1:duration_per_move/dt
            
            % %==============================================================
            % compute the local impulse response surrounding the center
            syn_sum=conv2(w_mask, SubMapLIP_fr); syn_sum=syn_sum(2:size(syn_sum,1)-1, 2:size(syn_sum,2)-1);
            Mask_lat_inh_eyepos=ones(3); Mask_lat_inh_eyepos(2,2)=0;
            Lat_inh_sum=conv2(Mask_lat_inh_eyepos, SubMapLIP_fr,'same');
            % Leaky integrator that simulation the dynamic change of the
            % local LIP_Eyepos neurons
            
            SubMapLIP= max(0,SubMapLIP + 5*(-SubMapLIP+syn_sum)*dt);                      %***
            act_record=[act_record max(max(SubMapLIP))];
            % Update the local activation
            LIP_Eyepos(-1+pos(1):pos(1)+1, -1+pos(2):pos(2)+1)=SubMapLIP;
            LIP_Eyepos=max(0,LIP_Eyepos);
            LIP_Eyepos_fr=mp2fr(LIP_Eyepos);
            if ~isequal(fig_opt, 'fig_off')
                subplot(5,2,7);mesh(LIP_Eyepos);ylabel('LIP Eye Position');  pause(pause_time);  %imagesc(LIP_Eyepos, clim_lip_ctr);
            end
            sim_clk=sim_clk+1;  
            %==============================================================
            
            % inhibit subitizing region, when target pos attended (note
            % inhibition taken place only when subitizing is still preceeding, it might disabled when no target found for subitizing ).
            PFC_NS_Input=zeros(1,size_PFC_N);
            max_nlip=max(max(NLIP));
            
% removied to fit scene integration            
%             % %++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
%             % Update only while subitizing
%             if  flag_subitizing==1  % in subitizing range
% 
%                 %DLPFC_NHiTuning=[0 .69 .83 .917 .96 .982 .987];
%                 PFC_NS_Tuning=[ .5 .6225 .7311 .8176 .8808 .9241 .9526];
%                 
%                 % start estimation while i>=3
%                 max_XNLIP=max(max(XNLIP));
%                 if i>=3    
%                     PFC_NS_InputNoise=random('norm', 0, noise_sigma, 1,size_PFC_N);
%                     PFC_NS_fr=exp(-((max_XNLIP-PFC_NS_Tuning+PFC_NS_InputNoise)/resp_sigma).^2);         %***
%                     % PFC_NS=max(0,PFC_NS+dt*(-PFC_NS + PFC_NS_Input-.5*SumLatInh));                     %***
%                     % PFC_NS_fr=mp2fr(PFC_NS);
%                     
%                     PFC_NR = max(0, PFC_NR+dt*(-PFC_NR +PFC_NS_fr-.5*SumLatInh));                                          %***
%                     PFC_NR_fr=mp2fr(PFC_NR);
%                     SumLatInh=PFC_NR_fr*(1-eye(size(PFC_NR,2)));
%                 end
%                 if ~isequal(fig_opt, 'fig_off')
%                     subplot(5,2,2);bar(num_count,PFC_NR);ylabel('PFC-N2');
%                     subplot(5,2,4);bar(num_count,PFC_NR);ylabel('PFC-N1');
%                     subplot(5,2,6); mesh(NLIP); zlabel('LIP-N');
%                     subplot(5,2,5); mesh(LIP_sal_fr);ylabel('LIP/FEF Saliency');pause(pause_time);
%                 end
% 
%                 %                     if ~isequal(fig_opt, 'fig_off')
%                             %                       subplot(5,2,6);mesh(NLIP);zlabel('LIP-N');
%                 %                 LIP_sal = max(0, LIP_sal + dt*5*(-LIP_sal + V2_out -
%                 %                 5*DLPFC_L2 - NLIP_inh));
%             %++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
            
            % %-------------------------------------------------------------
            % Update only while counting
%             else

                % %^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
                % Update activation of DLPFC spatial map
                % Layers DLPFC_L2 and DLPFC_L1 compose a recurrent network
                % which memorize the location of the items.
                % neurons in DLPFC_L2 is not easy to be activated, unless with the
                % cooperating act of LIP_Eyepos and hit(item recognized)
                % signal from IT.

                SynIn = LIP_Eyepos_fr+ceil(LIP_sal(center(1),center(2))-sal_hit_thresh);
                DLPFC_L1 = DLPFC_L1 + dt*(-DLPFC_L1 + 2*DLPFC_L2_fr+max(0,SynIn-DLPFC_Low_THRESHOLD));             %***
                DLPFC_L1_fr=mp2fr(DLPFC_L1);
                DLPFC_L2 = (DLPFC_L2 + dt*(-DLPFC_L2 + DLPFC_L1_fr));                      %***
                DLPFC_L2_fr=mp2fr(DLPFC_L2);

                if ~isequal(fig_opt, 'fig_off')
                   subplot(5,2,3); mesh(DLPFC_L1);ylabel('DLPFC L1'); %imagesc(DLPFC_L1, clim_dlpfc);
                   subplot(5,2,1);mesh(DLPFC_L2);ylabel('DLPFC L2'); %imagesc(DLPFC_L2, clim_dlpfc);
                end

                % Instantiate and update new layer in LIP, which responsible for cobine V2 output
                % and LIP_mem, and by which saccade direction is guided
                DLPFC_inh_mask=zeros(field_size);
                if (DLPFC_HitFromIT(1,1)==1)
                    msg='hit!';
                end
                % revised to increase inhibition range
                % DLPFC_inh_mask(field_size-pos(1)-mv+1,field_size-pos(2)-mv+1)=1;
                DLPFC_inh_mask(field_size-pos(1)-mv(1)-1:field_size-pos(1)-mv(1)+3,field_size-pos(2)-mv(2)-1:field_size-pos(2)-mv(2)+3)=1;
                DLPFC_inh=conv2(max(0,DLPFC_L2_fr-0.05), DLPFC_inh_mask, 'same');
                % %^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

% removied to fit scene integration                
%                 % %==========================================================
%                 % Update DLPFC_N
%                 % DLPFC_N_shift_floag 0: never shift; 1:begin to shift; 2:end shift
%                 if PFC_N_shift_flag==1
%                     PFC_NShift = conv(PFC_NR_fr, [0 0 1]);
%                     PFC_NShift = PFC_NShift(1,2:size(PFC_NShift,2)-1);
%                     shift_act=shift_amp*max(0, PFC_NShift);
%                     PFC_N_shift_flag=2;
%                 else
%                     shift_act=0;
%                 end
%                 SumLatInh=PFC_NC*(1-eye(size(PFC_NC,2)));
%                 PFC_NC=max(0,PFC_NC+dt*(-PFC_NC + PFC_NR_fr +shift_act));                  %***
%                 PFC_NC_fr=mp2fr(PFC_NC);
%                 PFC_NR=max(0,(PFC_NR+dt*(-PFC_NR +2*PFC_NC_fr  - SumLatInh)));                                          %***
%                 PFC_NR_fr=mp2fr(PFC_NR);
%                 if ~isequal(fig_opt, 'fig_off')
%                     subplot(5,2,2);bar(num_count,PFC_NC);ylabel('PFC-NC');
%                     subplot(5,2,4);bar(num_count,PFC_NR);ylabel('PFC-NR');pause(pause_time);
%                 end     
                %==========================================================
            %--------------------------------------------------------------
%             end
            
            % %++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
            % Updatae Common Activations
            LIP_sal =  max(0, LIP_sal + dt*15*(-LIP_sal + V2 - 3*DLPFC_inh));                    %***
            LIP_sal_fr= mp2fr(LIP_sal);
            
            if ~isequal(fig_opt, 'fig_off')
                subplot(5,2,5); mesh(LIP_sal_fr);ylabel('LIP/FEF Saliency');
                %subplot(5,2,3); imagesc(DLPFC_inh);pause(pause_time);
            end
        
            % Find the current max activation in the eye field
            % * has to be implemented as winner take all later
            pos_LIP=Extract2d_max(LIP_Eyepos);

% removied to fit scene integration
%             % Updata Numerosity map in LIP
%             NLIP= NLIP + 5*dt*(-NLIP +conv2(V2, NLIP_resp_mask, 'same'));                                     
%             XNLIP=1./(1+exp(-.5*NLIP));                                                                   %***
%             if ~isequal(fig_opt, 'fig_off')
%                 subplot(5,2,6); mesh(NLIP); zlabel('LIP-N'); 
%             end
            %++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

        end 
        
% %     Implement termination saccade mechanism
        tmp=zeros(field_size);
        tmp(pos,pos)=1;
        LIP_sal_tmp=max(0, (LIP_sal_fr - tmp));
        if(max(max(LIP_sal_tmp))<=sal_terminate_thresh)
            if stop_flag==1
                break; 
            else
                stop_flag=1; %stop in next run
            end
        end        
    end
    
    % Report number
%     [v idx]=max(PFC_NR);
%     n=min(6,idx-1);
%     t=i;
end

% Not used, a questionable funtion
function mv=SC(EyeField, center)
    pos=Extract2d_max(EyeField);
    mv =pos-center;
    mv(1)=mv(1)/max(1, abs(mv(1)));
    mv(2)=mv(2)/max(1, abs(mv(2)));    
end

% pp=Eye_field so skiped

function EyeField = update_eyefield(Img,pos)
    % create convolution vector for remapping
    field_size=size(Img,1);
    x=1+field_size-pos(1); % right ward saccade use left ward remapping filter
    y=1+field_size-pos(2);
    mask=zeros(field_size);
    mask(x,y)=1;
    
    % remapping
    EyeField = conv2(Img, mask); 
    l_bound=int16((field_size+1)/2);
    u_bound=field_size+int16((field_size-1)/2);
    EyeField = EyeField(l_bound:u_bound, l_bound:u_bound);
end


function pos=Extract2d_max(M)
    SizeM=size(M);
    Vec=reshape(M, SizeM(1)*SizeM(2), 1);
    [C, idx]=max(Vec);
    pos=[rem(idx(1), SizeM(1)),floor(idx(1)/SizeM(2))+1];
    pos=int16(pos);
end

function target_pos=Extract_NN(Act_map, center, threshold)
    Act_map=Act_map-threshold;
    [row col val]=find(Act_map>0);
    d_row=int16(row)-center(1);
    d_col=int16(col)-center(2);
    Dist=d_row.^2+d_col.^2;
    [v idx]=min(Dist);
    target_pos=[int16(row(idx)) int16(col(idx))];
end

function hold(delay_period)
start_time=clock();
end_time=clock();
    while delay_period>(end_time-start_time)
        end_time=clock();
    end
end

function fr=mp2fr(mp)
fr=1./(1+exp(-7*mp+3.5));
end